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1. Introduction 

Classical continuum theories lack any length scales and as such provide leading 
order approximations for a number of problems that contain microstructures. 
Microstructures typically introduce characteristic length scales that may induce gradient 
dependences of various kinds. Several continuum theories have been developed to 
deal with microstructures and their attendant phenomena. These include micropolar, 
micromorphic, strain-gradient, non-local, etc. See Eringen P and Nowacki [2j for 
catalogs of such phenomenological theories. Physical theories for complex fluids such as 
liquid crystals have also been introduced, based on symmetry breaking phase transitions 
that yield statistically-defined order parameters as additional thermodynamic variables. 
See de Gennes and Prost [3] for discussions of the fundamental principles of order 
parameter physics for liquid crystals. See Holm [3] for a variational description of 
order parameter theories of complex fluids. Often these phenomenological theories are 
combined with geometrical discussions based on the theory of Cosserat and Cosserat 

The mathematical theory of continuum mechanics for complex, or composite 
materials produced a number of interesting phenomenological models in the 1960's. 
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Among these are the models of differential type introduced by Rivlin and Ericksen jB] 
and the multipolar models of Green and Rivlin [7J [8J. The history of how these models 
were tested in comparison with experiments and refined by ab initio assumptions of 
thermodynamics is recounted, for example, in Dunn and Fosdick [§], Eringen jTj and 
Fosdick and Rajagopal [TO] . 

Recently, Bellout et al [IT| considered a fusion of the models due to Rivlin and 
Ericksen jH] and those of Green and Rivlin [71IH]. The resulting theory introduced higher 
order spatial velocity gradients into the energy that regularized the model solutions and 
endowed the model with promising stability characteristics. 

The present work specializes to a subclass of the Rivlin- Ericksen- Green multipolar 
fluids treated in Bellout et al ^JJ that has energy density given by 

\D\u\ 2 + DW(e, Ve, We, ...,D,VD, WD, . . .), (1) 

where D is the mass density and e = |(Vu + Vu T ) is the strain rate tensor. Materials 
whose energy density takes this form are called gradient fluids of degree n, where n is 
the order of the velocity gradients appearing in (^Q). The case n = is the Euler fluid 
(no velocity gradient dependence), while the case n = 1 coincides with the 2 nd grade 
fluid 0, whose energy depends on the velocity gradient through the strain-rate, e. 

Our aim here is to investigate the implications of adopting a subclass of these 
n th — gradient models for the well-known elliptic instability, which governs the rapid, 
violent transition from two-dimensional to three-dimensional motion at the onset of 
turbulence in Newtonian fluids [T2J E3 EH EES EH) • We shall not assess the implications 
for experimental measurements of this investigation, as we feel that such an assessment 
may still be premature. Instead, we continue the investigation begun by Bellout et al 
[TTj in studying the role of n th —gradient constitutive relations on fluid instability. We 
begin by casting the n th — gradient theory of nonlinear elasticity into the Euler- Poincare 
variational framework ^7j. The Euler- Poincare framework allows us to take advantage 
of several parallels between n th — gradient fluids and recently developed Lagrangian- 
averaged Navier-Stokes-alpha, or LANS— a, turbulence closure models of Foias et al [TH]. 
Since Rivlin JH] , remarkable parallels have been drawn between nonlinear elasticity and 
turbulence closure models. In our case, the Euler-Poincare framework leads to energy 
balance laws, a proper definition of momentum density, circulation theorems and to the 
Craik-Criminale (CC) class of exact solutions for the n th — gradient materials. 

The CC solutions [201 f° rm the basis for analyzing elliptic instability, in which two- 
dimensional flows with closed streamlines are subject to three-dimensional instabilities. 
Our aim in this paper is to determine the effects of n th — gradient viscoelasticity on the 
parametric resonance mechanism responsible for elliptic instability and on its growth 
rates. We follow the earlier treatment of elliptic instability for Newtonian fluids as 
reviewd, e.g., by Kerswell [TH], and we are guided by the results of Fabijonas and 
Holm |2~T| E2J based on the CC solutions for the LANS— a and similar closure models 
for turbulence. Thus, we consider plane wave disturbances of elliptical flows whose 
wave amplitude and wave vector are time-dependent. This approach leads to a Floquet 
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problem for the wave amplitude of the disturbance. Remarkably, we discover that these 

viscoelastic effects may be either stabilizing, or destabilizing, in the sense that they 

alter the shape and size of the instability domain while simultaneously increasing or 

decreasing the associated Lyapunov growth rates, depending on the parameter values. 

We hope that experimentalists may be guided by these results in testing whether 

n th — gradient models may be appropriate for the description of viscoelastic materials 

undergoing elliptic instability. 

The equations of motion for gradient fluids are obtained from the Eulerian form 

of Hamilton's principle introduced in Holm, Marsden and Ratiu called the Euler- 

Poincare theory for continua with advected quantities. In vector notation, this is 

d & L S L 

— m + u ■ Vm + (Vu) T ■ m + m div u — V — = , where m = — . (2) 

at 5D du 

For the class of Lagrangians we shall consider, for th , 1 st and 2 nd gradient fluids, one 

has 

L = J d 3 x{±D\u\ 2 + Du-R(x)-p(D-l) + LW(e, Ve)} . (3) 

The term in R in this Lagrangian boosts the gradient fluid flow into a frame rotating 
with angular frequency 2Q = curlR, while the term in p imposes the constraint D = 1. 
Hence V • u = 0, as implied by substituting D — 1 into the continuity equation, 

d t D + V ■ (Du) = . (4) 

Many mathematical regularity properties are available for the class of gradient fluids, 
especially for the case that the Lagrangian L in (JBJ) provides a norm (when evaluated on 
the constraint surface, D — 1). However, these regularity properties for gradient fluids 
will be discussed elsewhere, following Foias et al [2~3~j . 

The objective of the current paper is to investigate the stability properties of CC 
solutions of the gradient fluid equations. For CC solutions in an unbounded domain, 
the fluid velocity is linear in the spatial coordinate and the pressure is quadratic. The 
CC solutions may be regarded as the first term in a Taylor expansion in space, around 
a stagnation point of the gradient fluid flow in a moving frame. We shall use the theory 
of elliptic instability to investigate the exact nonlinear growth rates when CC solutions 
interact with a wave packet whose phase is frozen into the CC flow for gradient fluids 
of degree n = 1,2. (The Euler case n = was studied in the original work of Craik 
and Criminale See also Craik [21], Miyazaki [23] and Kerswell ^H] for subsequent 
developments. See also Lagnado and Simmen and Goddard and Alam (2Z| for similar 
analyses for an upper-convected Maxwell fluid and granular media, respectively.) 

Outline. Section|2]summarizes the properties of ideal gradient fluids that follow directly 
from their Euler-Poincare formulation. These properties include energy conservation, 
momentum balance and Kelvin circulation preservation, all of which follow from 
Noether's theorem. We then specialize to gradient fluids of degree n — 1,2. Sectional 
introduces the CC solutions for I s * and 2 nd gradient fluids. Section HJ discusses their 
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stability properties for both inviscid and viscous CC solutions. Here we introduce 
viscosity as in the theory of 2 nd grade fluids, to which the gradient fluids reduce when 
n = 1. Section El summarizes our conclusions. 



2. EP formulation of gradient fluids 

Hamilton's principle for first and second gradient fluids. The mathematical basis 
common to all ideal fluid motions is Hamilton's principle 

S J Ldt = 0, (5) 

where L is the Lagrangian for the system. We work in the Eulerian representation of 
fluids, where the Euler-Lagrange equation is replaced by the Euler-Poincare equation. 
See Ref. [T7] for a detailed discussion of Euler-Poincare theory. 

This paper focuses on the incompressible motion of first and second gradient fluids 
in a rotating frame. Thus, the class of Lagrangians we shall consider has the form 

L = J £(u, Vu, VVu, ...,D, VD, WD, R(x)) d 3 x , (6) 

where R(x) is the vector potential for the Coriolis parameter, i.e., curlR = 2f2. 
Specifically, we shall take, 

L = J d 3 x{±D\u\ 2 + Du-R(x) -p(D -1) +W(e,Ve)}. (7) 

Here p is pressure (a Lagrange multiplier), D is mass density, and u is fluid velocity. 
Through the function W(e, Ve), the first and second gradient fluids depend on e, the 
symmetric strain-rate tensor, 

L / Ou ■ du ■ \ 

§(«U + «i,i) = o 7U7 + 7d • (8) 



That is, first and second gradient fluids allow energy to depend upon strain-rate e 
and gradient of strain-rate Ve, respectively, PjJIHlE]]- The higher gradient fluids will 
allow energies that depend upon higher-order gradients of strain-rate. We introduced 
the dependence on e, Ve, etc., instead of Vu, VVu, etc., in equation (JJJ), so that 
the Lagrangian L will be invariant under rotations. Consequently, the resulting Euler- 
Poincare equations will admit an angular momentum balance relation and will satisfy 
the requirements of material frame indifference. 

Variational derivatives and natural boundary conditions. The variational derivatives of 
the Lagrangian J7J) for first and second gradient fluids are given by 

SL = J d 3 x{D(u + R) ■ 5u + (||u| 2 + u ■ R + W - p) 5D - (D - l)5p+ a : 5e|, (9) 

where a : 5e = tr(cr T - 8e) = a lJ 5eij, and we sum over repeated indices. The quantity 
a is the stress tensor, whose definition assures that it is symmetric, a T = a , 

5eij deij dSJeij ^ ^ 
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The added natural boundary condition for second gradient fluids, 
dW 

fi-^&„ = o, (ID 

arises from an integration by parts. Another application of integration by parts and use 
of the symmetry of gives 

J d 3 x{D{u + R - V ■ a) - 8u + (±|u| 2 + u ■ R + W - p)5D — (D — l)Spj, (12) 

where (V ■ cr)* = da 13 /dx 3 . Another natural boundary condition has been introduced 
and applied, 

n • a ■ 8u — , at the boundary. (13) 

This condition may be satisfied when the fluid velocity has no normal component at the 
boundary, by requiring that the normal stress have no tangential component, 

(n ■ cr) x n = , at the boundary. (14) 

The Euler-Poincare motion equation. The Euler-Poincare motion equation is 

d S L S L 

—m + u ■ Vm + ( Vu) ■ m + m div u — V-; — = , where m = — . (15) 
at oD on 

The momentum density m is defined as the variational derivative of the Lagrangian 

with respect to the fluid velocity u. For the gradient fluid Lagrangian ((7|), we see from 

(H2I) that this is 

8L 

m= — = D(u + R- V-cr). (16) 
ou 

We denote (Vu) T • m = nijVu J , and d/dt = d/dt + u ■ V is the material derivative 
along u. The incompressibility condition V • u = follows from the continuity equation 
d t D + V • (Du) = 0, evaluated for D = 1, as imposed by the pressure constraint. 
Consequently, the Euler-Poincare motion equation (fTo^) obtained from the gradient fluid 
Lagrangian is expressed as, 

j (u + R - V • a) + (Vuf ■ (u + R - V ■ cr) + V (p - \\u\ 2 - W - u • R) = , (17) 

together with V • u = 0. Next, we use the vector identity 

(u ■ V)R + (Vu) T ■ R = - u x curl R + V(u ■ R) , 

together with the Coriolis relation curlR = 2f2(x), and introduce the standard 
dissipation law for the first and second gradient fluids. Consequently, the motion 
equation takes the familiar form, 

— (u - V • cr) + (Vu) T • (u - V • cr) + 2fl x u + V (p - ||u| 2 - W) = uAu. (18) 
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2.1. Circulation theorem and energy-momentum conservation 

Kelvin- Noether circulation theorem. In the absence of forcing and dissipation, the 
Euler-Poincare theory for Lagrangians in the class (JOJ) provides a Kelvin-Noether 
circulation theorem J7j 

— i — — -dx = 0, 19 
at ,/ c ( u ) iJ ou 

which holds for integrations around any closed curve c(u) moving with the fluid. For 
the first and second gradient fluids considered here, this becomes, 

^- I (u - V ■ <x + R(x)) ■ dx = . (20) 
at J c(u) V / 

Stokes theorem then provides, for relative vorticity u) = curl u, that 

— / (w - curl(V • <t) + 20(x)) • = , (21) 
Js(t) v ' 

for any surface S(t) whose boundary dS(u) moves with the fluid. Consequently, we find 

the Helmholtz vortex dynamics equation for the total vorticity, in the absence of forcing 

and dissipation, as 

d 

— E + u- VE-E-Vu = 0, where E = u - curl(V • tr) + 2fi(x) . (22) 

Thus, the Kelvin-Noether circulation theorem in the Euler-Poincare framework implies 
that the total vorticity E is frozen into the flow of a non-Newtonian, first or second 
gradient fluid. Hence, its total vorticity E satisfies the Helmholtz vortex dynamics 
equation 

Energy conservation. From the Euler-Poincare theory, one may compute the 
Hamiltonian from the Lagrangian L in equation (J2J) for first and second gradient fluids 
in a rotating frame by applying the Legendre transformation, || 

H = (m,u) -L 

= J jd 3 x^D\u\ 2 + (T : e - DW(e,\7e) + p(D - 1)} - j h ■ a ■ udS (23) 

The corresponding conserved energy is found by evaluating this expression on the 
constraint manifold, D — 1, as 

E = J J d 3 x^\u\ 2 + (T :e-H/(e,Ve)| - j h ■ a ■ udS . (24) 

The surface integrals in the last two equations vanish, upon applying the normal-stress 
boundary condition (|14jl . for the situation in which the velocity u on the surface has 
no normal component. As a consequence, the inner product of the fluid velocity u with 
the motion equation (|TH|) yields, 

f \Vu\ 2 d 3 x + - (fh-V\u\ 2 dS. (25) 



dt J 2 . 

|| Actually, we compute only the Routhian; because we do not Legendre transform the pressure; and 
we do not complete the transformation to explicit dependence only on m. 
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The surface integral vanishes, in this energy balance relation for first and second gradient 
fluids in a rotating frame, provided u vanishes on the boundary. 

Momentum conservation. We express the Euler-Poincare equation ()15|) in components 



as 

d 5L d 



d 5L d ( 5L a SL _ • ^ SL n 

-u 3 ) + — -diU 3 - Ddi— = . (26) 



dt Su % dx- i \5u % / 5m J 5-D 
Observe that for a gradient fluid Lagrangian © given by 



£(u, Vu, VVu, . . . , D, VD, WD, . . .) d A x , 
we have variational derivatives 

d 2 dC 



SL 


dC 


d dC 


Id " 


ao 


dx 1 dD i 


SL 


dC 


d dC 


Su % 




dx 3 du % 4 



dx l dx m dD tm 
dx 3 dx 1 du % 



9 dC 

+ o / n i -+••• ( 27 ) 

where the — h • • • refer to any dependence of the Lagrangian density L on higher spatial 
derivatives of D and u. Therefore, upon performing the indicated differentiations by 
parts, one eventually finds the local conservation law for momentum, 

d $ L 

dtm-i = — — — tT] , with momentum density m, = -— , (28) 
ox 3 Su l 

and momentum-stress tensor T\ defined by 

T> = rn, W + (C- D^)SI - (« A + £- + ^ + " ' ') ■ <») 

Here we abbreviate, by using variational-derivative notation to denote, 

SL dC d dC d 2 dC 

+ n , n ^ ^r-u + ••• , 



6u k j du k j dx 1 du h jt dx l dx m du k - lm 
SL dC d dC d 2 dC 

■- + ■•-, (30) 



Su k jt du k 3l dx m du k jlm dx m dx n du k jlmn 
SL 8C d dC d 2 dL 



Suk jlm duk jlm dxnduk jlmn 9x n d X P d U k jlmnp 



+ 



The momentum conservation form (|28|) is guaranteed by the Euler-Poincare 
equation for any choice of Lagrangian that does not depend explicitly on the spatial 
coordinate. The Coriolis vector potential R(x) introduces explicit spatial dependence 
into the Lagrangian. Consequently, although not all components of the momentum will 
be conserved, we may still write the motion equation f)18|) as a momentum balance 
relation, 

d t mi = - ~^—T{ + z/Aui + e ijk u 3 2n k , (31) 

where the momentum-stress tensor T\ is given by ()29|) and is the completely 
antisymmetric tensor density, with e 12 3 = 1. 
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Equations (|27|) for the momentum density and (|29j) for the momentum- 
stress tensor indicate how the derivation and analysis may proceed within the 
Euler-Poincare framework for gradient fluids of degree three, four, five, etc. These 
generalizations correspond to allowing the strain-rate W in the Lagrangian J7J) to depend 
on higher gradients of the strain-rate e. Pursuing this direction further for n th degree 
gradient fluids is straightforward within the Euler-Poincare framework. However, the 
present paper stops at 2 nd degree gradient fluids. 

Choice of energy density W(e, Ve) for first and second gradient fluids. In this paper, 
we will examine elliptic instability via exact nonlinear Craik-Criminale (CC) solutions 
for specific cases that apply for first and second gradient fluids. For this study, we shall 
choose the strain-rate dependence in the potential energy density as a norm, 

W(e,Ve) = \a x \e\ 2 + \a 2 \S7e\ 2 , (32) 

where |e| 2 = e^e^ and |Ve| 2 = e^e^ in tensor notation. See also Bellout et al [1999] 
for a discussion of the role of this norm in proving the regularity properties of their 
Rivlin-Ericksen-Green multipolar fluids. The case ct 2 = 0, a± ^ corresponds to the 
equations for second gradient fluids, and ot\ = 0, a 2 = corresponds to the classic NS 
equations. For the choice in ()32|) . we have 

V • <j = aiAu — a2A 2 u. (33) 

Upon defining v = (1 — a\A + a 2 A 2 )u, the motion equation |TH|) takes the following 
form: 

d t \ + (u • V)v + (Vu) r ■ v + 2Q x u - z/Au 

+ V(p - ||u| 2 - |ai|e| 2 - \a 2 \ Ve| 2 ) = 0. (34) 

For this choice of the energy density, the stress tensor in (J29)) has the form: 

T? = niiUj +p5ij - (aie kj u k)i + ~a 2 (e jkt i + e k i d )u k ,h) 

= UiUj +p5ij - oil (eujUj + e kj u k)i ) + a 2 (eu Mm Uj - \{e jk + + e k ij)u k> u) . 

Note, this stress tensor is not symmetric. The Lagrangian in ((7|) is also not invariant 
under rotations, when the Coriolis vector potential R(x) is a fixed vector. In the absence 
of R(x), this Lagrangian regains invariance under rotations and the angular momentum 
in that case is conserved. However, the stress tensor in that case is still not symmetric. 



3. CC class of solutions for gradient fluids 

A solution to ()34j) on an unbounded domain may be obtained, by taking velocity in the 
linear form, u = S(t) -x + U(t) together with a pressure p , which is quadratic in space. 
The matrix S is a time dependent matrix such that 

S{j + Si m S m j + 2€i rnk H, m S k j = Mij, Sa = 0, (35) 
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and U(t) is the instantaneous velocity field at the origin. Here, M is a symmetric matrix 
defined as M^- = — <9j<9jP, where 

F =-J F-c/x + p (x,t) + (\J + S-V + 2Q x u) -x. (36) 

A typical approach is to choose a matrix S for which the left hand side of (|33j) is 
symmetric. Then, the corresponding pressure po( x , t) is determined a posteriori by (J36|) . 
We nondimensionalize the system using the variables x' = x//, £' = cut, u' = u/\u\l, 
v' = v/|o;|Z, a' x = af, = at 2 /l 2 , where / is a typical length scale and u = curlu . 

The resulting equation with the prime notation suppressed is (pHj) with v replaced by 
u/\u\. 

We construct a second solution to (JH4*|) of the form u + ui with corresponding 
pressure p + p x . We refer to u as the 'base' flow and U! as the 'disturbance.' The 
equations governing the disturbance are 

d t wi + u • Vvi + m ■ vi + (Vu ) T • vi + (Vui) T • vi + II x ui 

Ui • (V • o) - ||ui| 2 - ai(e )i 3 -(ei)ij - iai|ei| 2 



-V(pi 



«2(eo)ti,fc(ei)ij,fc - \ol 2 \ Vei| 2 ) = i/Aui, (37) 



with V ■ ui = 0, in which we mix tensor and vector notation, where v* = (1 — a\A + 
«2A 2 )uj and e,; = |(Vu,j + (Vuj) T ) for i — 0, 1, and II = 2f2 + curl v . In the above 
equation, we have used the fact that u is an exact solution to ()34j) together with the 
vector identity 

u ■ Vv + (Vu) T ■ v = curl (v) x u + V(u • v) (38) 
for any two vectors u, v. We choose the disturbance to be of the form 

ui = jLta(t) sin(/3 , 0(x, £)), (39) 
Pi = mi(t)cos((3ip(x,t)) +/i 2 pi 2 (t)cos 2 (/3^(x,t)), (40) 

^(x, t) = k(t) ■ x + S(t), and \i and (3 are scaling factors so that we can choose the initial 
conditions |a(0)| = 1 and |k(0)| = 1. The unknown phase ip(x,t) and the amplitudes 
a(t), pn(t), and pu^t) are to be determined. The incompressibility condition V ■ ui = 
gives 

a-k = 0. (41) 

Form this equation it follows that the nonlinear term Ui ■ v x in (|37j) vanishes exactly. 
Thus, in what follows, the sum uo + ui is an exact solution to the nonlinear equations of 
motion in (jMJ). By collecting on powers of sin(/3^>) and cos^VOj the evolution equations 
for the amplitudes and phase are 

p 12 - (T - l)|a| 2 + i/3 2 |a| 2 |k| 2 ( ai - a 2 \k\ 2 (3 2 ) = 0, (42) 

d t ip + (S ■ x + U) • k = 0, (43) 

d t (ra) +T5 T -a + n x a - (3Pk = -£jk| 2 a. (44) 
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Here , 

T(i) = 1 + «!/? 2 |k(t)| 2 + a#|k(t)| 4 , (45) 

-E'cj = vfi 2 is the vorticity based Ekman number, II = curlu + 2ft is the total 
vorticity of the system, and P = pn — ^/3 2 ai2L ■ (S + S T ) ■ k. Note that the amplitude 
scaling /x is immaterial. Without loss of generality, we set 

d t 5 + k • U = 0. (46) 

Then taking the gradient of (|43|) becomes 

d t k + S T ■ k = 0. (47) 

We eliminate the pressure term by taking the dot product of (jUj) with k and by using 
da/dt ■ k = —a • dk/dt = (S ■ a) • k, the first of which follows from (f4*Tj) and the second 
from (1471) : 

PP= J^jr[(S + S T ) -a] -k + n x a-k}. (48) 

In summary, we have obtained a new exact incompressible solution to (|34[). The variables 
are amplitude a(t) and wave vector k(£). Once these are determined, the pressure terms 
follow from (|42|) and (|4*Hj) . Note that u and u + Ui are exact solutions to the nonlinear 
equations, but Ui by itself is only a solution to ()34)1 linearized about uo. The exception 
is that in a rotating coordinate system (O ^ 0), ui is an exact solution by itself since 
this scenario corresponds to u = 1Z ■ x in a non-rotating frame, where 1Z is rigid body 
rotation about the z— axis; cf. Ref. |2HJ - We emphasize that the operator d t + S T acting 
on a vector represents the complete time derivative of that quantity in a Lagrangian 
frame moving with Uo- Finally, the equation for a(t) is^f 

d t 2L = y{2/3 2 (S ■ k) ■ k(«a + 2« 2 /? 2 |k| 2 )a - TS T - a- nxa + Pk- P w |k| 2 a|. (49) 
4. Elliptic instability for gradient fluids 

We examine the stability of a rotating column of fluid with elliptic streamlines whose 
foci lie on the ?/-axis: 

/ -1 + 7 0\ 
u = \uL • x, L— J 1 + 7 . (50) 

\ o oo/ 

Here, < 7 < 1 is the eccentricity of the ellipses, and the pressure is po = 
\oj 2 {1 — 7 2 )(x 2 + y 2 ). Equation (jTfjl with S = L is analytically solvable: 

k= [sin^costVl - 7 2 ), /tsin^sin^v 7 ! -7 2 ),cos#] T (51) 

where n 2 = (1 — 7)/(l + 7) and 9 is the polar angle that k makes with the axis of 
rotation. Equation ()49|) has the form 

d t a = Af(t; cci, a 2 , E m Q, 7, 9) ■ a, 

% Alternatively, one can collect on the terms linear and constant in x upon insertion of Uo + Ui into 
(|34|l . In either case, by enforcing 1)46(1 . both methods yield the same set of equations: ((42|l . (|47(l - (|49|l . 
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where the elements of the matrix M are periodic with period r = 2n / \J\ — 7 2 , the 
period of k(t). Therefore, the system can be analyzed numerically using Floquet theory 
[2H]- We compute the monodromy matrix V, that is, the fundamental solution matrix 
with identity initial condition evaluated at t = r. Equation (|44|) will have exponentially 
growing solutions if max, |3?(pi)| > 1, where pi, i — 1, 2, 3 are the eigenvalues of V, with 
corresponding Lyapunov-like growth rates given by 

a = ln{max |3?(pi)|}/r. 

i 

Thus, we can simulate numerically the solution to (|44ft over one period and indisputably 
determine the exponential growth rates. We can be certain that at least one of the 
eigenvalues will always be unity because of the incompressibility condition (}4~T|) and 
that the remaining two eigenvalues appear as complex conjugates on the unit circle or 
as real valued reciprocals of each other. 



4-1. Inviscid results for gradient fluids 

For flows with circular streamlines (7 = 0), the monodromy matrix can be analytically 
computed. It follows from (|51|) that |k(£)| = 1. Then, T is constant in time (denoted 
by T = 1 + ai/3 2 + a 2 /3 4 ) and (jUJ) has three linearly independent solutions: 

ai(t) = cos(£(t) + 0)k ±1 + sin(f (t) + 0)k ±2 (52) 

a 2 (t) = sin(£(t) + 0)k xl - cos(£(t) + </>)k ±2 (53) 

a 3 (t) = e z , (54) 

where = 2t(l + f2)cos#/Y , k_Q = [cos 9 cost, cos 9 sin t, — sin 8} T and kj_ 2 = 
[sint, — cost, 0] T are vectors orthogonal to k, and <fi is an arbitrary phase. Clearly the 
first two solutions a! and a 2 satisfy (}4"T]) . The monodromy matrix can be constructed 
from these three solutions: 

/ cos({(2tt)) cos#sin({(27r)) 0' 

V= -sin(£(27r))/cos# cos(^(2vr)) 

\ tan 0(1 - cos(£(2tt))) -sin^sin(^(27r)) 1 
The three eigenvalues are pi j2 = exp(±^(27r)), p 3 = 1. All of the eigenvalues lie on the 
unit circle, from which it follows that all solutions in the inviscid case for 7 = are stable. 
The values of cos 9 for which \p^\ — 1, i — 1, 2, 3 are called 'critically stable' and are given 
by £(27r) = run, m — 0, ±1, ±2, .... At these parameter values an exponentially growing 
solution can appear (together with an exponentially decaying one) as 7 increases from 
zero. Bayly jTHj argues that the evenness of Pk as a function of k implies that the 
eigenvalues, if real and unequal, must be positive. This dismisses the odd choices of m. 
Furthermore, Floquet theory is not applicable for the case m — 0. Thus, the possible 
choices for critical stability are £(27r) = 2nir, n = ±1, ±2, . . .. This corresponds to 

cos6= 2lTTny (55) 

These are the critical parameter value at which a(i) suffers exponential growth as 7 
increases from zero. For the NS equations (i.e. T = 1), only the n = 1 choice (called 
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the 'principle finger') is physically interesting. The other choices of n are extremely 
thin fingers with growth rates ten orders of magnitude smaller than that of the principle 
finger j22]. As ot\ and/or «2 increase from zero, however, the fingers widen and the 
associated growth rate increases. Finally, since |cos#| < 1, we conclude that there 
exists a band of stable eccentricities for 

-^<fi + l<^. (56) 

Additional understanding of this result emerges by following the analysis of Waleffe 
[To] and Kerswell ^H]- By taking the dot product of f)44j) with a, we obtain (for all 7 
and T) 

d (!N 2 ) „ 4 7 (T-1)M; 2| l2 

~^- L = -2T«i«2 + \ ; 7^yl*| 2 - (57) 

One can determine an exponential growth rate to leading order in 7 by inserting the 
zeroth order solutions for k and ai into the right hand side of this equation: 
1 d(jlaP) 
|a| 2 dt 

= - cos Of sin(2(£+ + 0)) - (1 + cos 6) 2 sin(2(£_ + 0)) 

- 2(1 - cos 2 6) sin(2t)] + - ^ sin 2 6 sin(2t), (58) 

where £± = £(t) ±t. Upon averaging over a period of a 1; this quantity will vanish except 
when £± = 0, corresponding to cos 9 = =fT /2(1 + Q). Compare this with (JHEJ). The 
maximum values for a will occur at <f> = =F7r/4 for £± = 0, respectively, with growth rate 

(2 + T ) 2 (2 + T + 2fi) 2 



- 7 + ( 7 2 ), (59) 



16 (2 + T ) 2 (l + fi) 2 

valid for To < 2 and Q not satisfying ()56jl . Thus, we see that the maximum growth rate 
increases as a function of oti and a 2 due to the T dependence of the critical stability 
point up to a maximum of a = 7, after which a set of stable solutions emerges in a band 
of nonzero eccentricities. See Fig. ^ 

For nonzero values of 7, we must investigate the system numerically. We use 
the variable coefficient ordinary differential equation solver DVODE jHOj- The level 
surface of the growth rate for fixed a.\ = is seen in Fig. El and Fig. 01 shows 
the growth rate surface maximized over the 7,cos6 l plane as a function of ai,Q!2- 
Numerical experiments show that cr max has the value associated with the NS equations 
for ai = a?2 = 0. As the parameters ai j2 increase, cr max increases to a value of unity on 
the line a\f3 2 + a^/? 4 = 1 + 2f2, and then decreases slowly to zero as ot\ t 2 — > 00. This 
threshold line corresponds to the maximal rate of change of cr max in ()59|) with respect 
to 7. See Fig. 3. 

4-2. Viscous results for gradient fluids 

The solutions to (|49J) must be simulated numerically for E u 7^ 0. An interesting feature 
of this equation is that, unlike the NS equations, a change of variables will not remove 
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T T 

Figure 1. Sample contour plots of instability regimes for a\ = 0.5, «2 = 0.3, j3 = 1.0, 
£7 W = computed on a 250 x 250 grid for various values of f2: (a) -0.5, (b) 0, (c) 1.0, 
and (d) 2.5. Note that the individual fingers touch the cosf axis according to 1)550. 

viscosity from the problem. However, the qualitative results for NS hold true here. 
Viscosity stabilizes the flow by lowering the maximum growth rate and introducing a 
stable band of eccentricities. This stabilization is slower than its NS counterpart, that 
is, the disspation in (J31j) is of the form uAu, not i^Av. See Figs. |U and El 

5. Conclusions 

The presence of nonlinear elasticity was seen to have profound effects on the properties 
of elliptic instability It can affect the growth rates, as well as the shapes and sizes of 
the unstable parameter regimes. One of the most profound effects is the thickening of 
the resonance domains (fingers, or Arnold tongues) in Fig. 1. These resonance domains 
of instability were predicted for the NS elliptic instability. However, in the NS case, 
they are infinitesimally thin. 

The second gradient fluid constituitive relation and the LANS— a turbulence model 
both introduce higher derivatives in the momentum density. We found that the highest 
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Figure 2. Surface of a = 0.01 for E u = 0, a% = 0, (3 = v2j = 0, and various 
a- 2 . Figure (a) shows the neutral surface for < ct2 < 0.25 and is an expansion of the 
boxed region in figure (b). For a.2 = 0, the critical stability point occurs at 9 = tt/3, 
which agrees with the classical NS results. The critical stability point shifts towards 
cos 9 = 1 as Q2 increases according to cos 9 = (1 + «2/3 4 )/2. As «2/3 4 exceeds unity a 
stable band of rotating flows with nonzero eccentricities appears. The corresponding 
surfaces for ai = and various ot\ are qualitatively the same, albeit as a function of 
ai/3 2 . See |2U- 








Figure 3. Growth rate surface a maximized over the cos 9, 7 plane as a function of 
a-y and <x<i for E u = 0, (3 = 1, (a) = and (b) = 2.5. In figure (a), the growth 
rate for ct\ = ui = (corresponding to the classic NS case) is 0.36. We see that the 
growth rate quickly increases to unity on the line ai(3 2 + a2/3 4 = 1 + 2f2 and then 
slowly decays. Not shown is that cr — > as a±, a.i — > 00. 

derivative dominates and produces qualitatively similar effects on the neutral stability 
surfaces. That is, Fig. 2 shows a similar behavior of the neutral surface as a function of 
CK2/3 4 as found for the LANS— a model as a function of «i/? 2 , in the present notation. 

As seen in Fig. 3, first and second gradient fluids increase the Lyapunov growth 
rates associated with elliptic instability for a\f3 2 + a^ 4 < 1 + 2f2 and then decrease the 
growth rates for parameter values beyond this threshold. When cti = 0, this relation 
recovers the result for LANS— a. Thus, the higher-order smoothing due to cti 7^ comes 
into play to reduce the maximum growth rate for short waves. 
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Figure 4. Surface of <x = 0.01 for = 0, (3 = 1, = 0, as functions of ai for (a) 
S w = 0.1, (b) E u = 0.5, and (c) = 1.0, JS^ = vjf* /\uj\. Again, the corresponding 
surfaces for a.\ = and various 02 is qualitatively the same. Namely, the unstable 
region shrinks as viscosity in E u increases. 




Figure 5. Surface of a = 0.01 for = 0, v/\u>\ = 1 as a function of the wavenumbcr 
(3 for (a) a2 = 0, a,\ = 1 and (b) oi\ = 0, a? = 1. For /3 -C 1, the flow is inviscid. As /? 
increases, the leading order term in an asymptotic cxapnsion of 1)44(1 will be — _E w |k| 2 a. 
Since E^ = vf3 2 /\Lu\, viscous dissipation quickly takes over. 

Viscosity has the expected effects on the domain of elliptic instability, as seen in 
Figs. 4 and 5. However, these effects depend sensitively on the value of ot\ and a 2 . 
Figure 5 shows how the effects of nonlinear mscoelasticity depend on the values of 
ai and c*2 as a function of the wave number f3. The 02 term corresponds to the /3 4 
dependence, which comes into play very rapidly in its effect on the neutral surface for 
elliptic instability in Fig. 5b. 

Our investigation followed the approach of Fabijonas and Holm [^TJ |22]> w ho 
studied the corresponding mean effects of turbulence on elliptic instability for a class 
of turbulence closure models. For inviscid fluids, the effects of elliptic instability seen 
in gradient fluids and in the turbulence closure models are qualitatively similar. The 
inviscid first gradient fluid corresponds to the LANS— a turbulence model, which can be 
viewed as the nonlinear terms in an LES model for turbulence whose filter is the inverse 
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of the Helmholtz operator (1 — «iV 2 ) ^H]- The inviscid second gradient fluid can be 
viewed similarly, for which the filter is (1 — aiV 2 + o^V 4 ) -1 , instead. 

Future studies may investigate the roles of other aspects of nonlinear stress on 
elliptic instability, for example, in the Rivlin-Ericksen-Green multipolar fluids analyzed 
in Bellout et al [TTj . 



Acknowledgments 

BRF thanks the Los Alamos National Laboratory for its hospitality, and gratefully 
acknowledges the financial support of the Laboratory's Turbulence Working Group and 
Theoretical Division. DDH is grateful for partial support by Imperial College London 
and by US DOE, under contract W-7405-ENG-36 for Los Alamos National Laboratory, 
and Office of Science ASCAR/AMS/MICS. 



References 



[1] Eringen A C, ed 1976 Continuum physics: Volume IV — Polar and Nonlocal Field Theories (New 

York: Academic Press) 

[2] Nowacki W 1986 Theory of Asymmetric Elasticity (Oxford: Pergamon Press) 

[3] de Gennes P G and Prost J 1993 The Physics of Liquid Crystals 2nd ed (Oxford: Oxford University 

Press) 

[4] Holm D D 2002 Euler-Poincare dynamics of perfect complex fluids, in: Newton P, Holmes P and 

Weinstein A, eds Geometry, Mechanics, and Dynamics (New York: Springer) 113-68 

[5] Cosserat E and Cosserat F 1909 Theorie des corps deformable (Paris: Hermann) 

[6] Rivlin R S and Erickson J L 1955 J. Rat. Mech. Anal. 4 323-425 

[7] Green A E and Rivlin R S 1964 Arch. Rat. Mech. Anal. 16 325-53 

[8] Green A E and Rivlin R S 1964 Arch. Rat. Mech. Anal. 17 113-47 

3] Dunn J E and Fosdick R L 1974 Arch. Rat. Mech. Anal. 56 191-252 

3] Fosdick R L and Rajagopal K R 1980 Proc. R. Soc. Lond. A 339 351-77 

1] Bellout H, Necas J and Rajagopal K R 1999 Int. J. Eng. Sci. 37 75-96 

2] Pierrehumbcrt R T 1986 Phys. Rev. Lett. 57 2157-9 

3] Bayly B J 1986 Phys. Rev. Lett. 57 2160-3 

4] Landman M J and Saffman P G 1987 Phys. Fluids 30 2339-42 

5] Waleffe F 1990 Phys. Fluids A 2 76-80 

3] Kerswell R R 2002 Annu. Rev. Fluid Mech. 34 83-113 

7] Holm D D, Marsden J E and Ratiu T S 1998 Adv. Math. 137 1-81 

3] Foias C, Holm D D, and Titi E S 2001 Physica 152-3D 505-19 

3] Rivlin R S 1957 Q. Appl. Math. 15 212 

3] Craik A D D and Criminale W O 1986 Proc. R. Soc. London A 406 13-26 

1] Fabijonas B R and Holm D D 2003 Phys. Rev. Lett. 90 124501 

2] Fabijonas B R and Holm D D 2004 Phys. Fluids 16 853-66 

3] Foias C, Holm D D and Titi E S 2002 J. Dyn. Diff. Eqns. 14 1-35 

i] Craik ADD 1989 J. Fluid Mech. 198 275-92 

5] Miyazaki T 1993 Phys. Fluids 5 2702-09 

3] Lagnado R R and Simmen J A 1993 J. Non-Newtonian Fluid Mech. 40 29-44 

7] Goddard J D and Alam M 1999 Part. Sci. Tech. 17 69-96 

3] Lifschitz A, Miyazaki T and Fabijonas B R 1998 Eur. J. Mech. B/Fluids 17 605-13 



EP formulation of n— gradient fluids 



17 



[29] Yakubovich V A and Starzhinskii V M 1967 Linear Differential Equations with Periodic 

Coefficients (New York: Wiley) 
[30] Brown P, Byrne G and Hindmarsh A 1989 SI AM J. Sci. Stat. Comput. 10 1038-51 



